/*!  \file  BdsQR.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      06/30/2013 11:07:54 AM
 *   \copyright GNU Public License.
 */
#include "BdsQR.h"
#include "Matrix.h"
#include "setThreads.h"
extern "C" {
void dbdsqr_(
const char   *  uplo,// 'U' or 'L'  B is upper/lower  bidiagonal
const int    *     n,// The order of the matrix B.  N >= 0.
const int    *  ncvt,// The number of columns of the matrix VT. NCVT >= 0.
const int    *   nru,// The number of rows of the matrix U. NRU >= 0.
const int    *   ncc,// The number of columns of the matrix C. NCC >= 0.
      double *     d,// On entry, the n diagonal elements of the 
							// bidiagonal matrix B.
         				// On exit , the singular values of B in decreasing order.
      double *     e,// On entry, the N-1 offdiagonal elements of the B
      double *    vt,// On entry, an N-by-NCVT matrix VT.
               		// On exit, VT is overwritten by P**T * VT.
const int    *  ldvt,// The leading dimension of the array VT.
      double *     u,// On entry, an NRU-by-N matrix U.
         				// On exit, U is overwritten by U * Q.
         				// Not referenced if NRU==0
const int    *   ldu,// The leading dimension of U.  LDU >= max(1,NRU)
      double *     c,// On entry, an N-by-NCC matrix C.
         				// On exit, C is overwritten by Q^T*C. 
         				// Not referenced if NCC==0
const int    *   ldc,// The leading dimension of the array C.
      double *  work,// dimension (4*N)
      int    * info);/* 
         = 0:  successful exit
         < 0:  If INFO = -i, the i-th argument had an illegal value
         > 0:
            if NCVT = NRU = NCC = 0,
               = 1, a split was marked by a positive value in E
               = 2, current block of Z not diagonalized after 30*N
                    iterations (in inner while loop)
               = 3, termination criterion of outer while loop not met 
                    (program created more than N unreduced blocks)
            else NCVT = NRU = NCC = 0,
                  the algorithm did not converge; D and E contain the
                  elements of a bidiagonal matrix which is orthogonally
                  similar to the input matrix B;  if INFO = i, i
                  elements of E have not converged to zero.
*/
}
BdsQR::BdsQR(): BiSVD("Lapack dbdsqr") { 
	aUpper    = true;
	aLower    = false;
	aUColMajor = true; 
	aVColMajor = true; 
	canBeParallel = true;
	setCompileTime(__DATE__, __TIME__);
}
bool BdsQR::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
void BdsQR::compute(const Matrix & B, Matrix & U, Matrix & Sigma, Matrix & VT){
	reset();
	assert( U.isColMajor()); assert(VT.isColMajor());
	
	if ( !B.isBiDiagonal() || !B.isSquare()) return;

	char uplo = B.isUpper() ? 'U' : 'L';
	int n = B.cols(), ldu = U.lda(), ldvt = VT.lda(), info = 0, lwork = -1;
	lwork = 4 * n;
	double * d = Sigma.getModData();
	double * e = NULL, * work = NULL; allocate(e, n); allocate(work, lwork);

	memcpy(d, B.getData()       , n * sizeof(*d));
	memcpy(e, B.getOffDiagonal(), (n-1)* sizeof(*d));
	U.makeIdentity();	VT.makeIdentity();

	int zero = 0, ldc = 1;
	dbdsqr_( &uplo, &n, &n, &n, &zero, d, e, VT.getModData(), &ldvt,
		 U.getModData(), &ldu, NULL, &ldc, work, &info);
	success &= info == 0;
	delocate(work); delocate(e);
}

REGISTER(BdsQR);
